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PACS 82 . 20 . -w - Chemical kinetics and dynamics 

PACS 03 . 65 . Vf - Phases: geometric; dynamic or topological 

PACS 05. 10. Gg - Stochastic analysis methods (Fokker-Planck, Langevin, etc.) 

PACS 05 . 40 . Ca - Noise 

PACS 87. 15. Ya - Biological and medical physics: Fluctuations 

Abstract. - We study a classical two-state stochastic system in a sea of substrates and products 
(absorbing states), which can be interpreted as a single Michaelis-Menten catalyzing enzyme or as 
a channel on a cell surface. We introduce a novel general method and use it to derive the expression 
for the full counting statistics of transitions among the absorbing states. For the evolution of the 
system under a periodic perturbation of the kinetic rates, the latter contains a term with a purely 
geometrical (the Berry phase) interpretation. This term gives rise to a pump current between the 
absorbing states, which is due entirely to the stochastic nature of the system. We calculate the first 
two cumulants of this current, and we argue that it is observable experimentally. 



Single molecule experiments [1] have led to a resurgence of work on the stochastic version 
of the classical Michaelis-Menten (MM) enzymatic mechanism that describes the enzyme- 
driven catalytic conversion of a substrate into a product [2]. A lot of progress (analytical 
and numerical) has been made under various assumptions and for various internal enzyme 
structures [3-7]. Still, many questions linger. Specifically, linear signal transduction proper- 
ties of biochemical networks, such as frequency dependent gains, are of a high interest [8-10]. 
These are studied by probing responses to periodic perturbations in the input signals (in our 
case, kinetic rates) due to fluctuations in chemical concentrations, in temperature, or due 
to other signals. However, such approaches usually disregard a phenomena known from the 
theory of stochastic ratchets [11], especially in the context of biological transport [12-15]. 
Namely, a system's symmetry may break, and, under an influence of a periodic or random 
zero-mean perturbation, the system may respond with a finite flux in a preferred direction. 
This, indeed, happens for the MM reaction [16], and we study the phenomenon in detail in 
this work. 

The ratchet or pump effect manifests itself during periodic driving of simple classical 
systems, such as a channel in a cell membrane, which is formally equivalent to the MM 
enzyme [16-26]. The prevalence of the transport terminology over the chemical one in the 
literature is grounded in experiments, where driving is achieved by application of periodic 
electric field or nonequilibrium noise that modulate barrier heights for different channel 
conformations. The resulting cross-membrane flux has contributions that have no analog in 
stationary conditions. 
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Fig. 1: Transition rates into and out of the absorbing states L (substrate) and R (product). Note 
that the rates fci and fc_2 denote total rates (rather than per mole of substrate/product), and 
they can be modulated, for example, by changing concentrations of the substrate and the product, 
respectively. 

In a related system, a quantum pump [27-29], pump fluxes admit a purely geometrical 
Berry phase interpretation [30]. A similar interpretation has been introduced for some 
special classical cases [31,32]. However, the Berry phase has not yet been derived for classical 
stochastic Markov chains that model chemical kinetics, even though recent developments 
strongly hint at the possibility [16,33-37]. For example, such systems admit introduction 
of a vector potential for the fluxes to characterize circular motion [34] , allow reformulation 
of the Langevin dynamics in the form that strongly resembles wave packet equations in 
quantum mecanics with a nontrivial Berry curvature of Bloch bands [36-38], and result 
in pump currents that depend on contour integrals over the the evolution of the system's 
parameters [16]. Although currents in most analyzed models were produced by the lack of the 
detailed balance rather than by external time-dependent perturbations, a dual description 
in terms of an external noise induced ratchet effect is usually possible. Conversely, time- 
dependent perturbations can break the detailed balance and induce the catalytic cycle [18]. 

In the present work we demonstrate that a theory based on the Berry phase can be 
constructed for a purely classical adiabatically slowly driven stochastic dynamics. The 
theory leads to an elegant interpretation of prior results, and it makes new predictions even 
for the minimal model of the MM reaction. Additionally, the theory fills in an important 
gap in the current state of the field by providing a simple, yet general recipe for calculation 
of mean particle fluxes and their fluctuations. 

In a standard minimal model, a single MM enzyme catalytic reaction consists of an 
enzyme-substrate complex formation and its decay into the enzyme and the product. Both 
reactions may be reversible [3]. We assume no multiple internal enzyme states (see [1,5] for 
the complementing discussion). The reaction can be viewed as a single bin (the enzyme) 
that can contain either zero or one substrate particles in it. The particle can escape from 
the bin by jumping into one of the two absorbing states, either returning to the substrate 
- the Left, or converting into the product - the Right, cf. Fig. 1. In turn, if the bin is 
empty, either of the absorbing states can emit a new particle into it. All transition times 
are exponentially distributed with the rates defined as in Fig. 1. We investigate the mean 
particle current J from L into R on time scales much larger than the typical transition times 
between any two states. Our main goal is to understand the effect of periodic driving, i. e., 
periodic rate changes. 

We emphasize again that this system is equivalent to a cross-membrane transport prob- 
lem, where the L/R states correspond to compartments on different sides of the membrane, 
and the bin is a membrane channel complex. Additionally, this system is relevant to trans- 
port through a quantum dot in the Coulomb blockade regime [39] . 

Let P e and Pf be probabilities of an empty/filled bin (unbound/bound enzyme). Ne- 
glecting particle discreteness and assuming that all rates change slowly, so that P e /f are 
always in equilibrium, one derives for the instantaneous L — > R current 



and the summation is over m = —2,-1,1,2. This is the classical (MM) current. Notice 
that, even if (kik^) = (k-ik-2) (where (. . .) denotes time averaging), J c \ = (j c i) ^ due to 
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nonlinearities of the MM reaction. 

How is this result affected by the particle discreteness? The master equation for the 
system's evolution is: 



(2) 



and we assume the adiabatic approximation, i. c., the rates k m oscillate with the frequency 
lo <C min m /c m . Since P e = 1 — Pf, (2) is equivalent to a formally solvable first order 
linear differential equation. However, expressing the solution in a simple form requires 
approximations (e.g., small fluctuations) even for a harmonic dynamics of the parameters 
[16]. Additionally, the solution provides little information about the fluxes beyond their 
mean values, and it is not easy to generalize to the case of more complicated, multi-state 
stochastic systems. Thus here we pursue a different analysis technique. 
The formal solution of (2) is 



P(t) = r (e- pfo). 
where T stands for the time-ordering operator, p(i) = [P e (t), Pf(t)] T , and 

H 



ki + k- 2 

— k\ - fc_2 
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(3) 



(4) 



It is useful to separate the parts of H that are directly responsible for the transition into 
or out of the i?-state, namely H = Hq — V + — V-, where 
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(5) 
(6) 



With this, we can write a formal expression for the probability P n to have n net transitions 
from the bin into R during time t. This is similar to [6], but we also allow fc_ 2 ^ 0, counting 
transitions from R into the bin with the negative sign. For example, 



P n =i = 1 

rt 



+ 



[ (ttiE/b(Mi)fcH(ti)E/b(*i,to) 

J to 
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Here 1 is the vector with all unit entries and 
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(7) 



(8) 



We introduce the probability generating function (pgf ) for the number of transitions from 
the bin into R, 



Z(X) 



(9) 



X is called the counting field, S(x) is the full counting statistics, and its derivatives give 
cumulants (or connected correlation functions) of P n , e.g. (n) = —idS{x)/dx\ x =Qi (S 2 n) = 
(n 2 ) - (n) 2 = (-i) 2 d2S(x)/d X \=o, etc. 
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We derive the full counting statistics similarly to [6,39], but we extend the method by 
allowing for an explicit time-dependence of H. First notice that (9) with the expansion 
of probabilities as in (7) is equivalent to an expansion of an evolution operator with a 
X-depcndcnt Hamiltonian: 

Z=l+f(e-^o 6{x ' t)dt ^p(t ), (10) 

where 



H( X ,t) = H (t) - V + (ty* - V_(t)e- ix 

fci + k- 2 -k-i - k 2 e lx 
-fei - k- 2 e~ lx fc_i + k 2 



(11) 



We define the two instantaneous eigenstates of H(x,t) as |«o) and \m). There are also 
the left eigenstates (u and (ui\ with the same eigenvalues; we normalize them so that 
(u n (t)\u m (t)) = 5 nm . 

Let's introduce an intermediate time scale St, 1/u) 3> St ^> max m (l/fc m ). During this 
time, typically, many transitions happen but k m 's are approximately constant. Consider 
time-points tj = to + jSt, j = 1 ... N, and i/v = t. In the adiabatic limit, one can ap- 
proximately rewrite the expression for the time ordered exponent in (10) as a product of 
evolution operators over time intervals of the size St with parameters set to constants at 
times tj, i. e., 

Z « l + e -H{x,t N )5t e -H(x,t N -i)5t e -#(x,to)« p ( to ). (12) 

Now we insert the resolution of the identity 

i = \uo{ti))(uo(ti)\ + |ui(ti))(ui(ti)| (13) 

after every exponent at i,. We define \uo(U)} as the eigenstate of H(x,ti) corresponding 
to the eigenvalue Ao with the smaller real part. Since St 3> 1/fci, | exp [— Ao(Xj U)St] \ 3> 
| exp [— Ai(x, U)5t] |, and terms containing \u\(t)) can be neglected. Moreover, after long 
evolution, any information about the initial and the final states is lost. Thus we rewrite the 
pgf as 

N 
i=l 

Finally, we approximate (u {ti)\uo{ti-\)) « cxp[— St(uo(ti)\dt\u (ti^i))}, which allows us to 
write 



T 



To 



T 



S(x)~S gcom + S cl = -— I dt[(uo\d t \u )+Xo(x,t)}, (15) 



where To = 2tt/oj is the period of the rate oscillations and T » To is the total measurement 
time. The last term in (15) would be the same even for a time-independent problem, and it 
has been discussed before [39]. Diagonalizing H, and denoting e± x = e ±lx — 1, we get for 
this term 



Scl= 2T~oJ o ° dt K ~ V R2 + 4 ( K + e +x + K - 



c) 



(16) 



The first term in (15) is more interesting. Since it depends only on the choice of the 
contour in the k space, but not on the rate of motion along the contour (at least in the 
adiabatic approximation), it has a geometric interpretation. We write 



<Sgcom — ,- 



r / A-dk, A ro = <u (k)|d fem K(k)>, (17) 
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where the integral is over the contour in the 4-dimentional parameter space (the k-space) 
drawn during one cycle. 

To estimate the error in our results, we note that (10) is exact, and then assumptions 
followed. First, we neglected the initial nonequilibrium relaxation on a time scale of ~ k^ 1 . 
Since the S(x) ~ T, this introduces an error in S of ~ l/(k m T), which vanishes for long 
observation times. Second, we projected the evolution only to the subspacc of states with 
the smallest real part of the eigenvalue. The resulting error is exponentially suppressed in 
the adiabatic limit by exp(— (Ai — X )/lu). Third, there is the coarse graining in (12). To the 
lowest order, it introduces errors in S in the form of commutators, such as [H(Xi £»)> ^(x> 
\ti-tj\ < St. Since (u (t + St)\[H( X , U), H( X , tj)]\u (t)} ~ 0{[ujSt} 2 ) for t < t u tj < t + St, 
this error is less significant than the 0(u>) contribution from the geometric term in S. Finally, 
the error due to the approximation of (uo(ij)|uo(ti-i)) in (14, 15) is of the same order. 

In a chemical system, the rates k m can be changed by many means, e.g., by varying 
the system's temperature. However, the simplest scenario is to couple the substrate and the 
product to particle baths and to vary the corresponding chemical potentials for both species. 
Since the rates k\ and fc_2 are proportional to the particle numbers, they will oscillate as 
well. Thus in what follows we assume that k\ and k-2 are time dependent, while ki and 
k-i arc constants. Then, using Stokes theorem, we write 

<[>A-dk= ff dfcirffc_ 2 F felife _ 2 , (18) 
where the integration is over the surface s c enclosed by the contour c, and 

We will call i 7 fe ll fe_ 2 the Berry curvature by analogy with similar definitions in quantum 
mechanics. The advantage of working with F rather than with the potentials A m is that 
the Berry curvature is gauge invariant, i.e., it does not depend on an arbitrary k-dependent 
normalization of |u (k)}. It is a truly measurable quantity. In our case, 

F - e- x (e ix fc2 + fe-i) 

kl ' k - 2 ~ [4«+e +x + 4«_e_ x + iP]3/2 • ^ 

Note that, with \ = 0, the Berry curvature is zero, so it is the counting field \ that 
introduces a nontrivial topology in the phase space of the eigenstates of H(x,t). More 
generally, a normal Markovian evolution corresponds to \ — 0, where the normalization of 
P ensures that F\ x=0 = 0. This may be the reason why nobody discussed the Berry phases 
in the context of Markov chains. 

Knowing the full counting statistics, one can study both the particle flux and its fluctu- 
ations. Importantly, from (15, 17), all cumulants will have the geometric and the classical 
terms. In particular, for the mean L — > R flux per unit time, we get 

where j c \ is the classical current defined in (1). Generally, the ratio of geometric and classical 
terms is ~ u>/k m <C 1. 

For the large observation time T, the total expected L — > R particle flux is (n(T)} = JT. 
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Similarly, its fluctuations are given by (S 2 n(T)) = J&T, where 
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T K 3 T K 5 
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K K 3 



(23) 
(24) 



The ratio of the pump and the classical noises is again ~ to/k m . Importantly, we see that, in 
general, expectations of the total, the classical, and the pump currents are not equal to their 
variances. Thus none of the currents in the MM problem is Poissonian. This is a general 
property of complex, multi-step reactions, which is often neglected in computational studies 

(21 

of biochemical reaction networks. Also note that J pump in (23) is not necessarily positive 
and can, in fact, decrease the total noise. However, the total noise variance J^ 2 ' is strictly 
positive. 

While the pump current in this system has been analyzed, cf. [16], to our knowledge, 
an expression for J parameterized by the rates, rather than by internal enzyme parameters 
has not been available. Even more importantly, expressions for the noise (either classical or 
pump) for variable kinetic rates have not existed either. 

The smallness of the pump effect compared to the classical current complicates its ob- 
servation. However, several opportunities exist. One is the dependence of J pU mp on the 
frequency of the perturbation, while J c \ ^ J c \(lo). The second possibility is to vary the rates 
along a contour with J c \ = 0. However, even in this case the noise may still be dominated 
by the classical contribution. 

To test our predictions, we choose one particular such contour with J c \ = 0, J pU mp ^ 
for a numerical analysis: 

ki = A + Pcos(wt), k- 2 = A + Rsm(ut), fc_i = k 2 = 1. (25) 

To estimate numerically, we implemented a Gillespie-like scheme [40] that admits an 

explicit time dependence of the rates. However, since J pU mp -C y J^i ■> suc h Monte-Carlo 
estimation of J would take too long. Instead, we numerically solved the master equation (2) 
for Pf(t) with a time discretization dt <C min m (l/fc m ). This gave a better precision than 
an analytical result of [16], which assumes R — ► 0. Knowing P/, j(t) = k^Pf — k-i{\ — Pf). 
Fig. 2 shows an excellent agreement of our theory with the numerical results for both J and 

J( 2 ). 

Intuitively, the flux cannot keep increasing linearly as 10 grows and the adiabatic ap- 
proximation fails. To understand the behavior at large frequencies, we consider the case of 
maxfc m <^ u). We look for the probabilities in the form P e /f = P e /f + SP e /f, where P e /f 
are calculated for the time-averaged parameter values, namely, P e = ((fc 2 ) + (k-i))/{K), 
Pf = 1 — P e , and SP e /f are small and oscillate fast. The latter can be found from (2). For 
example, when k\ = (k\) + 5fci cos(wi), fc_ 2 = (k-2) + <5fc- 2 sin(wt), and fc 2 , are con- 
stants, we find 5P e — —(SkiP e /ui) sin(a;t) + (<5fc_ 2 P//u;) cos(wi). Then the current averaged 
over the oscillation period is J = J c i({{k m )}) + Jpnmp, where J pump = <5fci<5fc_ 2 P e /(2tj) is 
the correction due to the fast rate oscillations. For the parameters as in (25), we again get 
J c i = 0, but J pump 7^ 0. However, the dependency of J pU m P on u) has changed from ~ ui in 
the adiabatic case to ~ 1/uj for fast oscillations. This agrees with [16], which used a different 
approximation and showed that there is a single maximum in the pump current at u) ~ k m . 

The phenomenon of </ pU m P 7^ can be explained in simple terms. Since a particle spends 
a finite time bound to the enzyme, the values of k\ and fc_ 2 cannot influence the system 
during the mean unbinding time following a binding event. If, during an oscillation, the left 
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Fig. 2: Comparison of analytical predictions and numerical results for a system defined by (25) 
with A = 1.5. Here the theory gives: J cl = 0, J pump = cji? 2 [2(2 - R 2 + 4A + 2 A 
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0. Points with error bars mark numerical results, while 



l-(50-2i? 2 )(25-2i? 2 ) 

lines correspond to theoretical predictions, (a) Mean L R particle flux J obtained by integrating 
(2) with dt = 0.001 and averaged over ~ 100 rate oscillations. Error bars are negligible, and so is 
the computational time to generate the data, (b) Mean L — > R flux variance for uj = 0.03 Hz. 
Each point is a result of simulating 3 • 10 6 elementary particle hops, and covers about 6000 rate 
oscillations. Error bars denote one standard deviation, as estimated from the posterior variance of 
j' 2 ' . These simulations took less than a day on a modern PC using Octave. 



binding rate is higher than the right one during the upramp of the cycle, then k\ "shields" 
growing values of fc_2 from having an effect, while fc_2 shields decreasing values of k\ during 
the downramp. This leads to a phase-dependent asymmetry that is the source of the pump 
flux. Larger frequencies lead to more shielding, hence the linear w dependence. However, 
for uj —> oo, the information about the phase of the oscillations is lost while a particle is 
bound, decreasing the asymmetry and the flux. 

Is the pump flux observable experimentally in biochemical (rather than channel trans- 
port) experiments? Exact zeroing of J c \, as in the numerical example above, which would 
make J pU mp the leading effect, may be difficult to achieve. However, the classical cur- 
rent is also small near the classical steady state, {k)\k2 = fc_i(fc_2). In this case, if 
fci,_2 = (ki-2) + <5fci,_2 cos(wt — 01,-2), and if the oscillations are small, 8k\-2/ki-2 <C 1, 
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then we can disregard variation of the Berry curvature inside the contour, and 

k 2 + k-i 



-(5fci<5fc_ 2 sin</>, (26) 



4um P ~ ^ -^^T L 5k i Sk -2 sin <j>, (27) 

where <f> = 4>- 2 — 4>i- This should be compared to the classical contributions in the same 
limit 

fc 2 <5fc^ — k-i8k 2 _ 2 + (k 2 — £;_i)(5fci<5fc_2 coscj) 
Jd « ^2 ( 28 ) 

(2) 2(fc 1 )fc 2 

We see that, while J pU mp/Jci ~ ^/(K), the ratio of the variances is further suppressed, 
^pum P /J cl ~ (u!/(K))SkiSk-2/(kik-2)- Thus overcoming the classical variance is the 
biggest concern for a successful experiment. 

Our model has a single substrate and a single product. Thus the enzyme we are discussing 
is an EC 5 enzyme. Properties of such enzymes vary dramatically depending on a reaction, 
a biological species, and mutations in protein sequences [41]. A typical range of k 2 in vitro is 
10~ 2 . . . 10 4 s _1 . Similarly, the Michaelis constant, Km = (fc-i + k 2 )[S]/ki in our notation, 
varies between 0.01 and 10 mM (here [S] is the substrate concentration). For our analysis, we 
take k 2 <~ 10 s _1 , and Km = 1 mM. While little is known about and k 2 separately, it's 
reasonable to assume fc_i ~ k 2 . Similarly, we take fc_ 2 ~ fci[P]/[5] since both rates are often 
dominated by the mean particle-enzyme collision time [3] (for example, for triose phosphate 
isomerase reaction, fc_ 2 /fci ~ 2[P]/[S'] [42]). Many enzymes with similar parameters have 
been characterized in the BRENDA database [41]. Then, with [5"] ~ [P] ~ 1 mM and with 

(5^1,-2/^1,-2 ~ 10%, we get Jp Um p ^ Jpump ~ 10" 2 /T , Jci - 10" 1 , and - 10 particles 
per second. Oscillation periods of ~ 1 s are attainable (and still satisfy u> <§C k m ), which 
puts the flux ratio at <~ 10 _1 . Different dependence on <f> can further improve detectability 
of Jpump compared to J c \. At these parameters, the pump flux becomes equal to the total 
flux standard deviation for an experiment lasting a few days. Alternatively, working with, 
say, ~ 10° enzymes, J pump and J pum p become observable in less than a minute. While 
real experiments will certainly have additional complications, it is clear that our predictions 
should be experimentally testable, at least in principle. 

In conclusion, we constructed the Berry phase theory of a purely classical adiabatic 
stochastic pump in a simple Michaelis-Menten enzymatic mechanism. Our approach allowed 
calculation of the particle flux and its variance, including the classical and the pump effect 
contributions. We believe that these predictions can be checked experimentally, and it should 
be interesting to consider their importance in the context of enzymatic signal processing. 

While we analyzed only one specific model system, the Berry phase approach is general 
and can be employed for many processes that can be reduced to slowly driven Markov 
dynamics. Examples range from charge transport in quantum dots at strong decoherence, 
to particle transport through cell membrane channels, and to various biochemical reaction 
systems. For these problems, our method has several advantages compared to already known 
techniques. First, it provides a formal recipe for an analysis of an arbitrary Markov chain: 
one should construct a Hamiltonian with counting fields, find its eigenvalue with the lowest 
real part and the Berry curvature from the corresponding eigenstate, and then the counting 
statistics is given by (15, 17). Second, the method provides a solution not only for average 
fluxes, but also for the full counting statistics and thus for all flux cumulants. This is 
important in the context of elimination of fast degrees of freedom for construction of coarse- 
grained biochemical reaction models (e. g., MM or Hill phenomenological laws), where a 
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correct treatment of noise in the remaining degrees of freedom has always been a point of 
contention. Third, the method allows to transfer many results of the well-developed Berry 
phase theory for dissipative quantum dynamics [43-45] to the field of chemical kinetics. In 
particular, techniques exist to extend our work and to calculate flux cumulants to all orders 
in u)/k min [46]. 
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